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Emergence and time evolution of micro-structured new-phase domains play a crucial role in determining 
the macroscopic physical and mechanical behaviors of iron under shock compression. Here, we investigate, 
through molecular dynamics simulations and theoretical modelings, shock-induced phase transition 
process of iron from body-centered-cubic (bec) to hexagonal-close-packed (hep) structure. We present a 
central-moment method and a rolling-ball algorithm to calculate and analyze the morphology and growth 
speed of the hep phase domains, and then propose a phase transition model to clarify our derived growth law 
of the phase domains. We also demonstrate that the new-phase evolution process undergoes three 
distinguished stages with different time scales of the hep phase fraction in the system. 

The high-pressure states of iron have long been of great interest because of its technological and sociological 
importance as well as its geophysical role within Earth core 1 . In particular, the structural phase transitions 
induced by pressure and temperature have been extensively studied 2 . The structural phase transition of iron 
under shock loading was first reported by Bancroft in 1956 3 , where the transition pressure was determined as 
about 13 GPa based on wave-profile analysis. After that, many research works have been done to investigate the 
transition pressure, transition mechanism and equation of state both in theory and at experiment 4 " 6 . In 2002, 
Kadau et at. first observed from large-scale molecular dynamics (MD) simulations the evolution process of phase 
transition from bec into hep structure in iron 7 . Later, they also studied the shocked polycrystalline iron 8 . In 2005, 
using the in situ X-ray diffraction technique with nanosecond resolution, Kalantar et at. directly confirmed the 
phase transition mechanism of iron 9 . That is, the bec-hep phase transition includes uniaxial collapse along the 
[001] direction and shuffling of alternate (110) planes of atoms. Using the same in situ technique combined with a 
modified Warren-Averbach method, in 2008, Hawreliak et at. derived the conclusion that single-crystal iron 
becomes nanocrystalline in shock transforming from bec to hep phase 10 , in reasonable agreement with results 
from large-scale MD simulations. 

In addition to its microscopic mechanism, clearly, a physical modeling of phase transition process for shocked 
iron also crucially requires a deepened knowledge about the nucleation rate, growth speed, and the associated 
morphology evolution of nanoscaled new-phase domains (here, hep domains), which keep unsolved up to now. 
This becomes particularly important when considering the above-mentioned experimental fact 10 that shock 
transformed iron is really characterized by nanoscale grains. Inspired by this observation, through systematic 
MD simulations and a rationalized theoretical analysis, here, we study the non-equilibrium phase transition 
process in iron under the critical pressure. For this purpose we develop a central-moment method and a rolling- 
ball algorithm to calculate and analyze the morphology and growth speed of the single hep phase domain. Then, 
we propose a phase transition model to help understanding our derived domain growth law. Finally, we dem- 
onstrate that the evolution process of hep phase follows a three- stage description with different time scales of the 
hep phase fraction in shocked iron. 

Results 

The sample material we use in simulations is single-crystal bec iron. The simulation tool is the well-known 
LAMMPS software package 11 . The interatomic interaction is described by an embedded atom method (EAM) 
potential 1213 , which has been proved to be able to successfully describe the mechanical properties and structural 
phase transition behaviors of iron under high pressure. The simulation box consists of 240 X 80 X 80 unit cells 
and contains approximately 3 X 10 6 atoms. The x, y, z axes are along the [100], [010] and [001] directions, 



SCIENTIFIC REPORTS | 4:3628 | DOI: 1 0.1 038/srep03628 



1 



respectively. The periodic boundary conditions are applied along the 
y and z directions to minimize the surface and edge effects. The shock 
wave compression is generated using the momentum mirror method 
along the x axis 14 . In order to study the non- equilibrium phase trans- 
ition process in detail, the shock velocities are chosen in between 300 
and 400 m/s, which produce pressures around the critical value for 
phase transition. The velocity interval between two successive simu- 
lations is 10 m/s. 

To identify the crystal structures, the atomic coordination num- 
bers and common neighbor analysis (CNA) values are calculated. 
The single phase domain atoms are extracted using the cluster iden- 
tification method. The phase interface shape is determined by our 
developed rolling-ball algorithm 15 . The growth speed of the phase 
domain is calculated using the central moment method and rolling- 
ball algorithm. 

Figure 1 shows the evolution snapshots from bcc to hep structure 
during the shock process. Figures l(a)-(e) show the phase transition 
mechanism. Here, the green and blue spheres denote the upper and 
lower layers of atoms, and the red and yellow spheres mean that the 
atomic displacements are along the [011] and [Oil] directions, 
respectively. After the shock wave swept samples along the [100] 
crystalline orientation, the atoms are compressed along shock dir- 
ection and form hexagons in the (011) and (Oil) planes, as shown in 
Fig. 1(b). Soon afterwards, some local atoms move through collective 
thermal fluctuations along [011] or [Oil] direction, as shown in 
Fig. 1(c). There is a relative slide between the two layers of atoms, 
which causes that the distance between each atom in the layers and its 
two second-neatest neighboring atoms along the y and z axes 
becomes farther. Once the defects form, they will drive the slip planes 
to slide alternately. The alternative slippage leads the bcc structure to 
evolve into the hep structure and the phase transition process is 
finished, as shown in Figs. 1(d) and 1(e). 

Figures 1(f) and 1(g) exhibit the formed hep phase domains in the 
shocked region, where the hep atoms are shown in blue color and the 
boundary atoms are shown in red color. From Fig. 1(f), it is clear that 
the nucleation sites of hep phase domains are randomly located. The 



initial morphology of the single phase domain is ellipsoid-like, see 
the inset to Fig. 1 (f) . There are only two kinds of phase domains in the 
system. One kind is to stack into ellipsoid-like along (011) crystal 
plane by alternative slippage along [011] and [Oil] directions. The 
other kind is to stack into ellipsoid-like along (011) crystal plane by 
alternative slippage along the [Oil] and [Oil] directions. With the 
formation and growth of the hep phase domains, different domains 
begin to interact and collide with each other. When the sliding types 
of two collided phase domains are the same, they link to form a bigger 
phase domain. If the sliding types of two collided phase domains are 
different, they form grain boundary and interact with each other, as 
shown in Fig. 1(g). 

In general, once a phase domain core forms, it will gradually grow 
up through the outward movement of the phase interface. The 
growth speed depends on the driving force of phase transition and 
the synergic movement of interface atoms, which are closely related 
to the shape of phase interface. Therefore, the morphology evolution, 
growth speed, and the interactions among phase domains have 
always been a focus in the study of phase transition kinetics. In this 
work, we extract the single phase domain atoms, directly determine 
and visualize the phase interface shape by detecting surface atoms of 
the phase domain, as shown in Fig. 2, where panels (a) -(c) show the 
atomic evolution of the single phase domain, and panels (d)-(f) are 
the corresponding evolution of the phase interface shape. One can 
observe from Fig. 2 that the phase domain forms ellipsoid-like by a 
superposition of several layers of atoms. The growth of the phase 
domain mainly has two ways. The first way is the slide of each layer of 
phase plane, as shown in Figs. 2(a) and 2(b). The interfacial atoms 
surrounding the phase domain change into hep arrangement via 
synergic movements. This type of growth is continuous. The second 
way is the successive activation process of new phase planes, as 
shown in Figs. 2(c), where two newly occurred phase planes are 
labeled as 13 and 14 in number. Compared to the slipping process 
of phase plane, the energy threshold for activating a new phase plane 
is higher, and an energy accumulation is needed. As a result, this type 
of growth is discontinuous. 
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Figure 1 | Microstructure evolution in iron from bcc to hep structure during the shock process. Panels (a)-(e) show the phase transition mechanism, 
while panels (f)-(g) show the formed new phase domains in the shocked region. 
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Figure 2 | Evolution snapshots of the single phase domain atoms and the corresponding phase interface shapes. Panels (a)-(c) are the atomic evolution 
of the single phase domain, while panels (d)-(f) are the corresponding phase interface shapes. Numbers in panels (a)-(c) denote layers of phase plane. 



In this paper, we calculate the central moment and principal- axis 
lengths according to the expressions (2) and (4) for multiple phase 
domains. The calculated results indicate that the three principal- axis 
directions (namely, three eigenvectors) approximately are [100], 
[Oil] and [Oil] for all phase domains. This demonstrates that the 
phase domain has different growth speeds along the shock loading 
direction, relative sliding and normal directions of phase planes. 
Figure 3 shows the principal-axis lengths and growth speeds of two 
phase domains which form at different times under the same shock 
velocity. One can observe that both the length and growth speed 
along the normal direction of phase plane are the smallest. The 
growth speed of the phase domain is prominently supersonic within 
a range 4 X 10 4 to 5 X 10 3 m/s. The time dependence of the principal- 
axis lengths can be approximately scaled as L ~ t 0A65 on average. In 
addition, we also note, from the growth speed curves along the [Oil] 
direction, remarkable oscillations with a period of —0.02 ps, which 
implies a discontinuous growth mode. The phenomenon is possibly 
caused by the fact that only the atomic layers exceeding a critical size 
are able to promote the activation of the new atomic layer. This is 
similar to the dislocation growth process, in which dislocation core 
should exceeds a critical size for the emission of dislocations. 

Furthermore, based on the above results, we calculated and plotted 
the time evolutions of surface areas and volumes of the two phase 
domains according to the surface area and volume formulae of ellips- 
oid, as the red symbols shown in Fig. 4. To verify the reliability of our 
proposed central-moment method, we counted the numbers of total 
atoms and surface atoms for each phase domain (blue symbols in 
Fig. 4), and calculated the surface areas and volumes of phase 
domains (green symbols) using our developed rolling-ball algorithm. 
From Fig. 4, it can be observed that the calculated results by the 
central-moment and rolling-ball methods are consistent. This fact 
also confirms that the shape of the single phase domain is always 
close to be ellipsoidal at early time. 

Theoretically, the evolution laws of surface area and volume of a 
single phase domain approximately follow A ~ t m and V ~ t n , 



Based on its ellipsoidal-like appearance, the single phase domain 
can be reasonably approximated as an ellipsoid with uniform density. 
Their central moments should be the same. By comparing the eigen- 
values and eigenvectors of central moments for the single phase 
domain and an ideal ellipsoid, we can determine the morphology 
parameters of the phase domain. The eigenvalues and eigenvectors of 
the central moment for the single phase domain satisfy 



Ax,- = fi-Zi (for i=l, 2, 3). 



(i) 



Here, the expression of the central moment A is given by 

n 

A=^m i (r,-f)(r i -r), (2) 



where m { and x { are the mass and position of the ith atom in the phase 

domain, respectively, and f = ^ J X{j n is the centroid. Whereas, the 

central moment of a uniform ellipsoid with the principal- axis lengths 
of a, b and c is given by 

/a 2 0 0 \ 
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where p is the atom density. Therefore, the three principal -axis 
lengths of the phase domain can be expressed as 
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Figure 3 | Principal-axis lengths and growth speeds of two hep phase domains which form at different times under the same shock velocity. Panels (a) and 
(b) are the principal-axis lengths, while panels (c) and (d) are the growth speeds. The symbols are the MD simulated results and the lines are the fitting results. 



respectively. If the principal- axis growth speed of the phase domain 
is constant, then m = 2 and n = 3. Otherwise, if the growth speed 
decreases with time, then m<2 and n < 3. The black curves in Fig. 4 
show our fitting results. We find that on average, the time evolution 
of the surface area and volume of a single phase domain can be 
respectively approximated as A ~ f 0 930 and V ~ t l 395 . The growth 
coefficients are almost invariant in the growth process of the single 
phase domain. But they are dependent of the nucleation times and 
positions of phase domains where the surrounding pressure are 
different. 

As a complementary clarification, now we propose a phase trans- 
ition model to illustrate the shocked kinetic process in iron based on 
the order parameter theory of Ginzburg- Landau. For this purpose we 
choose the slippage (0 of the lattice as the order parameter, which 
varies from £ = 0 in bec structure to £ = 1 in hep structure, as 
schematically shown in Fig. 5(a), where the horizontal axis represents 
the distance away from the phase interface. In general, for the uni- 
form (bulk) phase transition of iron, the system that undergoes 
transformation from bec to hep structure through the lattice slippage 
needs to overcome a potential barrier 16 , as schematically shown in 
Fig. 5(b). However, for the shocked iron, a solely uniform description 
is insufficient and the phase domain effects should be reasonably 
included in a phase transition model. In the nucleation stage of the 
phase domains, the atoms in the local region of stress concentration 
overcome the potential barrier by collective thermal fluctuations. 



From above numerical results, we have obtained that in the growth 
stage of a phase domain, the growth speed is supersonic and the stress 
wave has no time to propagate in the hep phase domain. Therefore, 
the growth of a phase domain is mainly driven by the interface 
energy. Figure 5(e) shows the potential energy distribution of a slice 
in the simulated system, where the regions of red, blue, and other 
colors represent the bec, hep, and interface structures, respectively. It 
is obvious that the potential energy in the interfacial region lies in 
amplitude between those in bec and hep regions. Thus, it is now clear 
that the interface energy is negative and prominently reduces the 
potential barrier between two phases, as schematically shown in 
Figs. 5(c) and 5(d). As a result, the transition process becomes easier. 

With keeping this physical picture in mind, the energy of the 
system expressed with order parameters reads approximately 



F = 



D , 



/(f) 



d 3 r, 



(5) 



where /(O and (V£ 



energy of the system, respectively. /(£) = -f 



are the bulk free energy and interface 

1 



■f is a 



3 

bi-stable function with two stable points, £ = 0 and { = 1. Here, the 
parameter a is a system parameter related to temperature and pres- 
sure. Under low pressure, a > 111, the bec structure is stable, while 
under high pressure, a < 1/2, the hep structure is stable. The possible 
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Figure 5 | The potential energy distribution and schematic diagram of phase transition. Panels (a)-(d) are the schematic diagram of phase transition, 
where the horizontal axis represents the distance away from the phase interface. Panel (e) shows the potential energy distribution of a slice in shocked iron. 
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anisotropy in the interface energy has been ignored for simplicity. 
The evolution equation of the order parameter can be expressed as 

d t z=-=m-DVt. (6) 

For the steady growth of one-dimensional phase domain, £ = £(rj) = 
£(x — c 0 t), and it satisfies the following eigenvalue equation 

f'(Z)-D%Z + cod n Z = 0, 
4—oo=0 ? (7) 

I ^7 — ^ + go ~ 1? 

where c 0 is the growth speed of the hep phase domain. To describe the 
growth of a three-dimensional phase domain, we adopt the local 
coordinate system instead of the Cartesian coordinate, r = r 0 + Xn 
+ fiti + vt 2 , where r 0 represents a point at the interface, and n, t l5 1 2 
represent the normal and two principal tangential unit vectors of the 
interface at the position r 0 , respectively. The evolution equation of 
order parameter can be rewritten as 

d t £ =/'({) ~D(d\ + dl + d 2 v + (h +k 2 )d x -k 1 dp- k 2 d v ) £, (8) 

where k Y and k 2 are the curvatures along the two principal tangential 
directions, respectively. According to the relation £ = £(rj) = ^{k — 
v£), the above expression can be reduced to 

f(0-Dd 2 ^+(-Dk+v)d^ = 0 (9) 

with the boundary conditions £ |^_oo = 0, £ |^ +00 = 1, and k = ky 
+ k 2 . Comparing Eqs. (7) and (9), the growth speed of phase domain 
can be evaluated as 

v = c 0 + Dk. (10) 

For shock-induced phase transition, the interface energy is related to 
the pressure surrounding the phase domain, and D is a function of 
pressure. From above simulated results and theoretical analysis, we 
have obtained that the growth speed is supersonic and the D is almost 
invariant in the phase domain growth process. Therefore, from Eq. 
(10) we get that the growth speed of the phase domain is a function of 
the local curvature. When the volume of phase domain is initially 
small, the local curvature is large and the energy that induces phase 
transition is relatively more concentrated. Thus, the growth speed is 
relatively high. With the growth of the phase domain, the interfacial 
area becomes much larger and the local curvature decreases, Therefore, 
the energy for phase transition becomes more dispersed and the 
growth speed decreases. This is consistent with our MD results. 

For an ellipsoidal phase domain, the local curvature is non- 
uniform on the surface of the phase domain, and the larger the 
curvature, the higher the growth speed. Therefore, the phase domain 



becomes more and more flat or prolate with time. Actually, it has 
been shown in Fig. 1(g) that various phase domains have evolved to 
be disc-like, spherical, columnar, elongated, etc., in the later stage of 
shock loading. For the moment it is helpful to give a very simple but 
illustrative estimate on growth dynamics of the phase domain. For 
this purpose, if the phase domain is regarded as a sphere with radius 
R, the growth speed Eq. (10) reduces to R = c 0 + 2D/R. When 
hypothesizing c 0 — > 0, one obtains that the linear length of phase 
domain is R(t) = V^Dt. Interestingly, our MD simulation results, 
which show that the linear length of an ellipsoidal phase domain is 
L — f A65 , close to this analytical expression. The difference in the 
exponents (0.465 versus 0.5) is caused by the difference in the cur- 
vatures between a sphere and an ellipsoid. 

The material properties, such as the constitutive relation and 
equation of state, are significantly influenced by the transition frac- 
tion of the new phase 17 . In the non- equilibrium phase transition 
process under shock, the mechanical behaviors are coupled with 
the phase transition process. The phase transition fraction is relevant 
not only to the evolution of a single domain but also to the interac- 
tions among neighboring domains. Therefore, the phase transition 
fraction is a highly concerned quantity in the studies of material 
phase transition. Figure 6 (filled circles) shows the MD results of 
the time evolution of phase transition fraction under two different 
shock velocities. For every shock velocity, remarkably, there are three 
obvious stages, which accordingly represent different evolution char- 
acteristics of phase transition fraction. By comparing with atomic 
evolution images of phase transition process, we find that in the stage 
between points A and B (the first stage), new phase domains succes- 
sively form and each phase domain independently grow up. In the 
stage between points B and C (the second stage), nearly no new phase 
domains form and the existing phase domains further grow up to 
interact with each other. In the stage after the point C (the third 
stage), the phase domains get saturated in the perpendicular direc- 
tions and grow up only along the shocking direction. 

Theoretically, the phase transition fraction is a function of nuc- 
leation and growth rates of the phase domains. The nucleation num- 
ber I follows approximately I — t m , while the atomic number AT of a 
single phase domain follows approximately N — t n . In the case that 
both the nucleation of new phase domains and the growth of the 
existing phase domains simultaneously happen in the system, then 
the phase transition fraction / follows approximately/ — t m+n . We 
have fitted our MD results, see Fig. 6. On average, we find that the 
evolution of phase transition fraction scales approximately/— t 1 89 in 
the first stage,/ — f 123 in the second stage, and/— f 080 in the third 
stage. 

Discussion 

Shock-induced phase transition process of iron from bec to hep 
structure has been investigated via systematic molecular dynamic 




time(ps) time(ps) 



Figure 6 | Phase transition fraction under two different shock velocities. Solid dots denote the MD results, while the red, green and blue lines represent 
the fitting results of three different evolution stages. The shock velocity is set at 380 m/s in panel (a) and 400 m/s in panel (b). 
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simulations and theoretical analysis. Based on the shape character- 
istics of the hep phase domains, we have proposed a central-moment 
method to calculate and analyze the morphology and growth speed 
of a single phase domain. It has been manifested that in the initial 
independent growth stage, the domain morphology is ellipsoid-like 
with three principal axes approximately along [100] , [01 1], and [Oil] 
directions. The growth speed of a single phase domain is prominently 
supersonic within a range 4 X 10 4 to 5 X 10 3 m/s. We have shown 
that on average the size, surface area, and volume of the single hep 
domain have their time evolutions L ~ f A6 \ A ~ f - 930 , V ~ t l 39 \ 
respectively. Based on the order parameter theory of Ginzburg- 
Landau, we have presented a phase transition model to explain our 
found growth law of the single phase domain. Finally, we have 
demonstrated a three-stage evolution law for the phase transition 
fraction in shocked iron. 

Methods 

Three methods have been employed in our numerical simulations and data analysis: 
(i) The numerical MD simulations are performed using the well-known LAMMPS 
software package. The interatomic interactions used in the simulations are described 
by embedded atom method potentials. The shock wave compression is generated 
using the momentum mirror method; (ii) The atoms are distinguished by the com- 
mon neighbor analysis (CNA) method. In this method the signature of the local 
crystal structure around a selected atom is identified by computing three character- 
istic numbers for each of the N neighbor bonds of the central atom. The single phase 
domain atoms are extracted using the cluster identification method; (iii) The phase 
interface shape is visualized by the rolling-ball algorithm. The growth speed of the 
phase domain is calculated using the central- moment method and rolling-ball 
algorithm. 
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